
fs = 1000;         % 采样频率 1kHz
N = 1000;
sigSourceCfg.freqs = [100, 200, 300]; % 复正弦信号的频率（Hz）
sigSourceCfg.gains = [1, 0.8, 0.6]; % 复正弦信号的频率（Hz）
sigSourceCfg.noisePower = 1e-2; % 噪声功率
sigSourceCfg.fs = fs; % 采样率（Hz）
sigSourceCfg.numSamples = N; % 样本数
% 生成信号
x = sinu_generator(sigSourceCfg);

p= 10;

% period
h_per = fftshift(fft(x,N))/N;
S_per = h_per.*conj(h_per);
S_per_db = pow2db(S_per);

% AR 谱估计
% [ar_coeffs, NoiseVariance] = aryule(x, p); % 使用 aryule 函数估计 AR 系数
% [ar_coeffs, NoiseVariance] = ar_yule(x, p); % 使用 ar_yule 函数估计 AR 系数
[ar_coeffs, NoiseVariance] = ar_cov(x, p); % 使用 ar_cov 函数估计 AR 系数
% 计算功率谱密度
freqs = linspace(-pi, pi, N); % 频率范围
H = freqz(sqrt(NoiseVariance), ar_coeffs, freqs); % 计算频率响应
psd = (abs(H).^2) / (2*pi); % 计算功率谱密度
psd_db = pow2db(psd);

% ideal
% H_ideal = freqz(sqrt(sig_pow), a, freqs); % 计算频率响应
% psd_ideal = (abs(H_ideal).^2) / (2*pi); % 计算功率谱密度
% psd_ideal_db = pow2db(psd_ideal);

% 绘制功率谱密度
fig = figure('Visible','off');
plot(freqs, S_per_db); hold on;
% plot(freqs, psd_ideal_db);
plot(freqs, psd_db);
legend('period','ar');
title('AR Spectrum Estimation');
xlabel('Frequency (rad/sample)');
ylabel('Power Spectral Density(dB)');
exportgraphics(fig,'test.jpg');
close(fig);